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Abstract 

We study a stochastic epidemic model consisting of elements (organisms in a community or 
cells in tissue) with fixed positions, in which damage or disease is transmitted by diffusing agents 
("signals") emitted by infected individuals. The signals decay as well as diffuse; since they are 
assumed to be produced in large numbers, the signal concentration is treated deterministically. 
The model, which includes four cellular states (susceptible, transformed, depleted, and removed), 
admits various interpretations: spread of an infection or infectious disease, or of damage in a 
tissue in which injured cells may themselves provoke further damage, and as a description of 
the so-called radiation- induced bystander effect, in which the signals are molecules capable of 
inducing cell damage and/or death in unirradiated cells. The model exhibits a continuous phase 
transition between spreading and nonspreading phases. We formulate two mean-field theory (MFT) 
descriptions of the model, one of which ignores correlations between the cellular state and the signal 
concentration, and another that treats such correlations in an approximate manner. Monte Carlo 
simulations of the spread of infection on the square lattice yield values for the critical exponents 
and the fractal dimension consistent with the dynamic percolation universality class. 

PACS numbers: 05.50.+q, 05.70.Ln, 87.10.Hk, 87.23. Cc 



* e-mail: |fernandopereira bh@gma il.com| 
'''e-mail: dickman@fisica.ufmg.br 



I. INTRODUCTION 



In many epidemic-like processes disease spreads via an agent emitted by the affected 
elements (cells or organisms) themselves [l| • Epidemics have been modeled extensively using 
deterministic reaction-diffusion equations [2j, and stochastic particle systems 00. In the 
simplest epidemic models, such as the susceptible-infected-removed (SIR) and susceptible- 
infected- removed-susceptible (SIRS) processes 0, disease is transmitted by contact between 
infected and healthy organisms, without explicit representation of a transmitting agent. But 
since the latter may have a dynamics of its own, typically involving diffusion and decay, it is 
of interest to include this agent explicitly in a more complete description, particularly when 
the spatial structure of the epidemic is analyzed. 

A similar observation applies to a viral infection, and to the spread of damage in tissue 
following irradiation. In the latter case, initially affected cells may become sources of a signal 
that damages nearby cells, which were not harmed in the initial event. These secondary 
cells may then become additional sources of the harmful signal. Such a scenario has been 
proposed for the radiation-induced bystander effect (RIBE) 8NlO| . While the direct result is 
damage to some or all of the irradiated cells, the long-term effect is characterized by damage 
or death of unirradiated cells or "bystanders". It is thought that irradiated cells release a 
signal (possibly a cytokine) that diffuses through the medium, causing damage in previously 
healthy, unirradiated cells Thus signal molecules in RIBE play a role analogous to a 
disease agent in an epidemic. 

In this work we introduce an epidemic model in which damaged or infected elements 
briefly emit signals; the latter diffuse and decay at a certain rates. Healthy organisms or 
cells may become infected or damaged due to the presence of the signal, and so may become 
new sources. We formulate the model on a discrete two-dimensional space, the square lattice. 



Epidemic models with spatia 



intensively in recent years 
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structure and short-range interactions have been studied 
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LH14|, and applied to the spread of disease in humans 
and plants 15l-ll7l|. Here we focus on processes initiated in a single cell or organism, or 
in a localized region. Key questions are then (1) the rate of spreading, as reflected in the 
growth in the number of affected individuals, and their spatial distribution, and (2) whether 
spreading continues indefinitely, limited only by the size of the available region, or stops 



before attaining a size comparable to the of the system. In the infinite system-size limit, 
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the two regimes are separated by a phase transition. In the supercritical (spreading) phase, 
there is a nonzero probability to spread indefinitely, whereas in the subcritical (nonspreading) 
phase the process dies out with probability one. 

Phase transitions in stochastic epidemic models with spatial structure have received con- 
siderable attention; an important example is the general epidemic process (GEP) 

el a. 

The GEP is essentially a stochastic susceptible-infected-removed (SIR) model with spatial 
structure. Initially, all individuals are susceptible (S) except for one or a few infecteds (I). 
Susceptibles with one or more I neighbors become infected at a certain rate A, while infecteds 
recover at rate /i, after which they are forever immune, hence removed (R) from interactions 
with other individuals. Transmission (S+I — > 21) is typically restricted to nearest-neighbor 
S-I pairs on a regular lattice or a network. The GEP exhibits a phase transition as the ratio 
X/fi is varied. The supercritical phase is characterized by a growing active region, which 
invades regions containing susceptibles and leaves behind an inactive region composed of 
individuals in states S and R. Activity is thus restricted to a "ring" separating as yet unaf- 
fected and formerly active regions. (In a finite system, the final state is completely inactive.) 
Analysis of the GEP shows that its critical behavior belongs to the dynamic percolation uni- 
versality class (3), Jjj]. If the process is modified so that a recovered individual can become 
susceptible (SIRS model), it is possible to maintain an active stationary state in which the 
processes of infection, recovery, and loss of immunity occur continuously. The SIRS model 
with spatial structure exhibits a phase transition belonging to the directed percolation (DP) 
universality class jsllisj], exemplified by the contact process 19, 20]. 

In this work we study a model in which individuals may be in one of four states: sus- 
ceptible (S), transformed (T), depleted (D), and removed (R), the latter class designating 
individuals that have died or are otherwise sequestered from the rest of the population. The 
principal new feature of our model is the mechanism by which infection is transmitted: the 
transition from susceptible to transformed is mediated by a signal released by cells in state 
T, rather than via direct contact. Such cells may recover (becoming susceptible once again), 
may be removed, or may emit signals. In the latter case, the cell immediately enters state 
D, after which it may either recover or be removed. Although cells in states T or D may 
recover, there is a finite probability of permanent removal. Thus we expect that, as in the 
GEP, the process will exhibit spreading and nonspreading phases, with activity concentrated 
in a ring. Assuming the phase transition is continuous, it is reasonable to expect the critical 
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behavior to be that of dynamic percolation. Given the novel mode of transmission, it is 
nonetheless of interest to verify this assumption. 

The remainder of this paper is organized as follows. We define the model in Sec. II, 
and in Sec. Ill discuss two mean-field approaches, a simple one that ignores diffusion, and 
a more detailed formulation that takes diffusion into account while still assuming spatial 
homogeneity. In Sec. IV we present simulation results for the phase diagram, critical 
behavior, cluster properties, and spreading velocity. A summary and a discussion of our 
results are provided in Sec. V. 

II. MODEL 

The model is defined on a square lattice of L 2 sites, each of which hosts an individual (an 
organism or a cell, depending on the choice of interpretation). Individuals exist in one of 
four states: susceptible (S), transformed (T), depleted (D), or removed (R). In addition to 
the discrete state variable, each site bears a signal concentration CV,- > 0. Individuals 
emit signals upon making the transition from state T to state D; we adopt concentration 
units such that each such event produces one unit of signal molecules. The transitions 
between states £11*6 cLS follows: 

(i) An individual in state S, at site has transition rates [iCy and uCij to states T 
and R, respectively. These are the only rates that depend on the signal concentration. 

(ii) An individual in state T has transition rates wts, wtd and wtr, to states S, D and 
R, respectively. 

(iii) An individual in state D has transition rates w^s an d wda, to states S and R, 
respectively. 

The states and transitions are summarized in Fig. [TJ note that state R allows no escape. 

The signal concentration Cy(t) evolves via diffusion and decay. We assume that the 
number of signal molecules is very large, so that the concentration may be treated deter- 
ministically, via the equation 
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T -D 

FIG. 1: States and allowed transitions. The rates for transitions out of state S are proportional to the local 
concentration of signal. The latter is produced when individuals move from state T to state D. 



— ^ = VACij - XCij + 

k=l 

where A denotes the discrete Laplacian, D is the diffusion rate, A is the decay rate, riy is 
the number of times site has made the transition from T to D, and the t^-ij are the 
transition times. Since the n^j and the transition times are random variables, the Cjj are 
as well. 

In the limit of very low signal concentration, the discrete nature of the signal molecules 
makes an important contribution to the fluctuations. Thus our continuum, deterministic de- 
scription may not be suitable for the low-concentration limit. Another possibly troublesome 
aspect of the diffusion equation is that, given a localized source at time zero, the concen- 
tration at any time t > is nonzero (albeit very small) at points arbitrarily far from the 
source. While this is unphysical, we do not expect it to cause any significant effect in the 
system under study. Indeed, the diffusion equation is widely applied, with apparent success, 
in modeling biological transport, and systems of reaction-diffusion equations (including ap- 
propriate noise terms) have been found to yield reliable predictions for critical properties of 



nonequilibrium systems 18 



The model involves a rather larger set of parameters: the coefficients /i and u, the diffusion 
and decay rates, and five additional transition rates. It is nevertheless clear that large values 
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of \i and wtd/(wtr + wts), an d small values of A, favor spreading. Evidently, spreading 
can only occur for T> > 0. Note that for wds — 0, there is no functional difference between 
states D and R, and we have a simpler, three-state model. 

We are primarily interested in an initial condition consisting of the origin in state D, with 
an associated unitary signal concentration, and with all other sites in state S, and free of any 
signal. The ensuing evolution corresponds to an epidemic spreading from the origin. The 
current size of an epidemic can be defined as the number of individuals in states other than 
S. Of particular interest are the number Nn(t) of removed individuals, the number iVr(i) of 
individuals in state T, and the (spatial) average signal concentration, C(t). The latter two 
quantities reflect the degree of spreading activity, since, if both are zero, further advance 
of the epidemic is impossible. In the spreading phase, starting from a small, localized set 
of affected individuals, Nr and C grow without limit in an unbounded system, whereas 
in the nonspreading phase these quantities cease to grow after a certain time. In a finite 
system, Nr and C, must eventually saturate, even in the spreading phase. The distinction 
between spreading and nonspreading phases is nonetheless evident in large, finite systems 
since, in the spreading phase, a finite fraction of individuals are eventually affected, whereas 



in the nonspreading phase the final fraction of affected individuals is ~ 1 /L 2 2l|, |22| . It is 
worth noting that, strictly speaking, an absorbing configuration corresponds to one without 
transformed cells, and with the signal concentration everywhere zero. Since C decays at a 
finite rate, such a situation can only obtain in the infinite-time limit. The implications for 
defining survival are discussed in Sec. IV. 

For simplicity, we assume that the signal-dependent transitions (i.e., from S to either T 
or R) have rates that are proportional to the local signal concentration. Other dependencies 
are conceivable; at the end of the next section we briefly consider rates proportional to C 2 . 



III. MEAN-FIELD ANALYSIS 



In the simplest mean-field analysis we factorize the joint probability distribution for the 
iV-site system into a product of single site probabilities, and treat the signal concentration 
Cij as independent of the state of the site. Denoting the probabilities for site (i, j) to be in 
state S, T, D, or R by Sy, Tij, Dij, and Rij respectively, and the mean signal concentration 
by C^, we then have: 
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dC 

- V (Cm - Qj) + wtdTh - XQj (2) 



dt 



dSij 

~dT 

dTij 



dt 
dD 



ij 



dt 
dRij 

dt 



(ki,ij) 

-(/! + rfCijSy + wrsTij + wosDij (3) 

fiCijSij - (w TS + w TD + w TR )Tij (4) 

wmTij - (w DS + w DF t)Dij (5) 

vCijSij + wmTij + w DR D i:j (6) 



where the sum in the first equation extends over the nearest neighbor sites (k, I) of site 
If we take the continuum limit of these equations and let X(r) = (C, S, T, D, R), we obtain 
a set of reaction-diffusion equations, dX/dt = DV 2 X + f(X), in which only the element D cc 
of the diffusion matrix is nonzero. 

We study a spatially uniform mean-field theory, which corresponds to the fast-diffusion 
limit. In this case the MF equations for the site probabilities are as above (removing the 
subscripts ij on all variables), while the concentration satisfies 

^ = wT-\C (7) 
dt 

Given the large set of parameters, it is convenient to fix all but one, which then plays 
the role of a control parameter. Somewhat arbitrarily, we choose w = wtd (i-e., the rate 
at which transformed cells emit signal and become depleted), as the control parameter, and 
denote its critical value by w c . We use the uniform analysis to set basic limits on survival 
of the spreading process, by studying an epidemic scenario. That is, we consider a very 
small initial source probability D(0), S(0) = 1 - D(0) ~ 1, and T(0) = R(0) = 0. (We 
set -D(O) = C(0) = 10~ 13 .) Then, as t — > oo, a non-vanishing value of R corresponds to an 
epidemic in which a nonzero fraction of individuals are affected, i.e., to the spreading phase. 

Figure [2] shows the evolution of the transformed fraction T(t); in the nonspreading phase 
T[t) decays monotonically, while in the spreading phase it grows at intermediate times. The 
depleted fraction D(t) exhibits a similar behavior. The fraction of removed individuals, 
R(t), grows steadily in the spreading regime, until it saturates; in the nonspreading regime 
it saturates at a very small value (see Fig. [3]). In the spreading phase, the growth regime is 
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followed by a crossover to exponential decay, marking the end of the epidemic. The crossover 
occurs due to the depletion of susceptibles. 

A simple analysis allows us to determine the boundary between spreading and nonspread- 
ing phases. Since we are interested in the early stage of the evolution (following the initial 
transient) we set S = 1, and seek a solution in which the probabilities T and D, and the 
signal concentration grow exponentially: T(t) = Tie 11 , and similarly for D and C . Inserting 
these expressions in the MF equations yields 

(7 + A) (7 + w T ) = fiw (8) 
where wt = w + wts + wtr- Equating 7 to zero yields the critical threshold: 

\{w T r + wts) , n s 



Note that spreading is impossible for [i < A, and that w c is independent of parameters v, 
Wos an d wdr- The above conclusions are verified in numerical solutions of the (uniform) 
MF equations. 




100000 200000 

t 



FIG. 2: Transformed fraction T(t) in uniform mean-field theory. Parameters A = 0.05, fi = 0.2, v = 0.1, 
wtr = wdr = wts = wjjs = 0.2, and (lower to upper) w = 0.13, 0.13333... = w c , and 0.134. 



The original MF equations, Eq. ([6]), not only treat sites as independent, but also treat the 
concentration at a site as independent of its state. This is a rather drastic approximation, 
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FIG. 3: Uniform mean-field theory: limiting fraction of removed individuals versus w. 



because signal molecules are only created when a site undergoes the transition T — > D; 
other sites only acquire a nonzero signal concentration via diffusion. This approximation 
can be improved by introducing a concentration variable Cj for each site type; here Cj 
denotes the mean concentration at a site, given that it is in state J. To derive a set of 
mean-field equations for the site probabilities and the associated concentrations, we treat 
the amount of signal at a given site as a discrete variable, n. Let P(J,n,t) denote the 
(joint) probability that a site is in state J and has exactly n units of signal. Then the 



probability of state J is P( J, t) = J2 n P{J-> n i Oj an d Cj(t) = nP(J, n, t)/ J2 n P{J-> n i t) = 
^ n nP(J,n,t)/P(J,t), so that 



Consider, for example, state S. There are contributions to dP(S,n,t)/dt due to (1) decay 
of the signal; (2) diffusion between the site and its neighbors; and (3) transitions between 
state S and other states. To treat diffusion at this level of approximation, we suppose that all 
four neighbors of a given site have the same, average concentration, C(t) = P(J, t)Cj(t). 
Then we have 



dC 



dt 



E n n[dP(J, n, t)/dt] [dP(J, t)/dt]Cj(t) 
P(J,t) P(J,t) 



(10) 
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dP ^ n ^ = \[(n + l)P(S,n + l,t)-nP(S,n,t)} 

+ ADC[P(S, n - 1, t) - P{S, n, t)] + AD\(n + 1)P(S, n + l,t)- nP(S, n, t)] 
+ w TS P(T, n, t) + w DS P(D, n, t) - n(fi + u)P(S, n, t). (11) 

Summing over n, we find, 



P(S,t) 



P(T) P(D) 
+ ^tsCtj^ + Wd8 c d ^> } (12) 



where {n 2 ) s = ^2 n n 2 P(S,n.t)/P(S,t). The second term in Eq. (|T0|) involves, 

Z n [dP(S,n,t)/dt] _ (a + u)c+w P (T) +W C P{D) 

— --(/. + v)C s + w TS — + WDsC D p^. (13) 



Multiplying by C$ and subtracting the result from Eq. (|T2l) . we obtain 

^ = _AC 5 + AV[C - C s ] + ^t S [Ct - <?s] 59 + ^s[C D - C^^y, (14) 

where we have set var(n) = (^ 2 )s — C 2 S to zero, as is usual in a mean-field approach. 
Proceeding in the same manner, we find, 

= _XC T + AV[C - C T ] + fiCs[C s - Qr]^§j, ( 15 ) 

= -\ Cd + AV[C -C D ]+w[l + C T - C D ]^r, (16) 

and 
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^ = -XC R + 4V[C-C R ] + uC s [Cs-C R }^ 

+ w TR [C T -C R }^ + w DR [C D - C R ]^. (17) 

Numerical solution of this set of equations yields a critical threshold, w c , which decreases 
monotonically with diffusion rate T>, approaching the value of the simple MF analysis as 
T> — > oo. Figure S] shows the evolution of the state probabilities and concentrations in a near- 
critical system, as predicted by the diffusive mean- field theory (DMFT). The predictions of 
DMFT for w c are compared with simulation results in the following section. 
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FIG. 4: Diffusive mean- field theory: temporal evolution. Lower panel (lower to upper): state probabilities 
P(D,t), P(T,t), total signal concentration C(t), and P(R,t); upper panel (lower to upper): conditional 
concentrations Cs and Ct (indistinguishable at this scale), Cr, and Cr>- Parameters A = 0.05, fi = 0.2, 
v = 0.1, wtr = Wdr = w ts = w ds = 0.2, T> = 0.02, and w = 0.2, slightly greater than w c = 0.19587. 



An advantage of MFT is that it readily allows investigation of diverse scenarios. We use 
MFT to perform a preliminary study of a variant of the model defined above, in which the 
rates for the transitions S — > T and S — > R are \iC 2 and uC 2 , respectively, representing a 
situation in which healthy individuals are essentially insensitive to very low signal concentra- 
tions. In the epidemic context, this corresponds to a situation in which small concentrations 
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of a disease agent are effectively eliminated by the immune system, whereas higher concen- 
trations overwhelm it. Experience with nonequilibrium phase transitions in systems such 



as Schlogl's second model 



23|, leads one to expect a discontinuous phase transition in this 



case. This is indeed verified for certain sets of parameters, featuring relatively large values 
of n; an example is shown in Fig. |HJ (Note that in this case the transition value w c depends 
on the initial signal concentration.) 
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FIG. 5: MFT: Final removal fraction Roo versus w for fi = 10, C(0) = 0.01, when the transition rates from 
S to T and S to R are oc C 2 , with other parameters as in Fig. [2] jumps from about 0.16 to 0.998 when 
w = 0.15335. 



IV. SIMULATIONS 

Simulations of the model defined in Sec. II are performed on square lattices of L x L 
sites. These studies begin with all sites in state S (and having signal concentration zero) 
except for the origin, which is in state D, with a signal concentration of unity. Although 
the boundaries are open, the system size is chosen such that sites on the boundary remain 
in state S, and have negligible signal concentration, for the duration of the studies. In the 
simulation algorithm, at each step, corresponding to a time interval At, the cellular states 
evolve as follows: 
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(i) if site is in state S, it remains in that state with probability 
p s = exp[— (fi + u)CijAt}. It makes a transition to state T with probability (1 — p s )fi/ (fi + u), 
and to R with probability (1 — p s )u/(fi + v). 

(ii) if the site is in state T, it remains in this state with probability p t = exp[— wxAt]. 
The probability of a transition to state J (=S, D, or R) is (1 — pt)wrj 

(hi) if the site is in state D, it remains so with probability pd = exp[— woAt]. The 
probability of a transition to state J (=S or R) is (1 — Pd)wuj /wd- 

(iv) sites in state R remain in this state. 

In addition, the signal concentration is updated in accord with Eq. (CQ), that is, CV, — > C'^ = 
exp(—\At)Cij + UAtJ2ki[Cki — Cij], with the additional contribution CV,- —> C^ + l at a step 
in which site makes a transition from state T to state D. Most of the studies reported 
below use At = 0.4. We found that the results using a smaller time step (At = 0.2) we 
essentially the same. In studies of large V values, however, it was necessary to reduce the 
time step to 0.2 or 0.1 to maintain reliability. 

Since the transition between spreading and nonspreading phases is apparently continuous, 
we determine the critical point w c by searching for power-law behavior of the level of activity 
n(t), the survival probability P(t), and the mean-square distance of removed cells cells from 
the origin, R 2 (t), as functions of time. The activity level n(t) is conveniently defined as the 
number of sites in state T, since growth of the active region requires transformed cells or 
a nonzero signal concentration. We find that at long times, the behavior of the total signal 
concentration is similar to that of Nt- 

As noted in Sec. II, the definition of survival is not completely obvious in the present 
model. Our definition is based on the continued increase in the number of R cells. To 
begin, we study the distribution of waiting times t w between successive events in which 
N R — > Nr + 1. Studying a large number of realizations up to some maximum time, t^, we 
find that they can be divided into two classes: those in which Nr increases until the end, and 
those in which this number saturates well before. In the latter class, the final configuration 
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is devoid of T cells, and the total signal concentration is extremely small, so that further 
growth in Nr is essentially impossible. We find that in the class of realizations with sustained 
activity, the probability distribution of waiting times, P(t w ), falls to zero above a certain 
value. On this basis, we define a maximum waiting time time, t^ M , somewhat larger 
than the value associated with the cutoff of P(t w ). If, in a given realization, the waiting 
time t w exceeds twM, the system is taken to be inactive, and the realization ends; the 
time of deactivation is taken as that of the most recent transition to state R. The survival 
probability P(t) is then defined as the fraction of realizations that are active at time t. (For 
the parameter set analyzed below, we use t\v m — 180.) 



The quantities P(t), n(t), and R 2 (t) are expected to follow 
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n(t) ~ t\ (18) 
Pit) ~ t~\ (19) 
R 2 (t) ~ t Zsp , (20) 

at the critical point, where 9, 5 and z sp are critical exponents associated with spreading. 
We expect the cluster generated by the critical process to be a fractal; its radius of gyration 
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should follow R g ~ n 1 / ? , where Df is the fractal dimension 

We perform detailed studies using the parameter values A = 0.05, /i = 0.5, v = 0.1, and 
wtr = wdr = wts = wds = 0.2, using system sizes L of up to 600, and simulation times of 
up to tu = 10 000 time units. (The number of simulation steps is t^/At.) For each value of 
the diffusion rate D studied, we determine the critical value w c using the power-law criterion 
for P(t). Specifically, we estimated w c using the local slope 0(t), as described below. The 
uncertainty in w c , on the order of 5 x 10~ 4 , reflects the range of w values for which we cannot 
rule out an asymptotic linear behavior of 6(t) versus 1/t at long times. The resulting phase 
boundary is compared against mean-field theory in Fig. El Mean-field theory furnishes the 
correct order of magnitude, but underestimates w c , especially for small values of T>. The 
diffusive MFT furnishes a slight improvement over simple MFT. For T> — > oo, the mean-field 
prediction converges to 2/45 = 0.0444...; simulations in this limit (using a spatially uniform 
signal concentration) yield w c = 0.045(1), consistent with MFT. Figure [7] shows a similar 
comparison for w c as a function of /j,, for T> = 0.02; in this case the DMFT prediction is 
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about a factor of five smaller than the simulation value. 

We also determined w c for a rather different set of parameters: v = 0.51, \x = 0.79, 
A = 0.1, V = 2.94, and w^r = w^s — w tr = wts = 0.0079. In this case simulation 
yields w c = 3.15(5) x 10~ 3 , while simple and diffusive MFT yield w c = 2.29 x 10~ 3 and 
2.31 x 10~ 3 , respectively. The closer agreement between simulation and MFT in this case 
can be attributed to the higher diffusion rate. 
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FIG. 6: (Color online) Critical rate w c versus diffusion rate V for A = 0.05, [i = 0.5, v — 0.1, and 
wtr = Wdr = wts = wds = 0.2. Points: simulation; solid curve: diffusive mean-field theory. The dashed 
line indicates the value predicted by simple MFT. Error bars are smaller than the symbols. 



A. Critical behavior 

Following preliminary studies which indicated that w c ~ 0.297, (for the parameter set 
mentioned above, with T> = 0.02), we performed a more detailed study using L = 360 
and tu = 5 000, with 72 000 realizations for each value of w studied. A large number of 
realizations is necessary to obtain a clear result for w c and the critical exponents. Using a set 
of 10 3 or 10 4 realizations (a number that would be sufficient for studying the contact process, 
for example) the results are dominated by fluctuations. We believe that this is due to the 
multi-step nature of propagation, and in particular, to diffusion. For the relatively small 
diffusion rate used here, the initial stages of propagation depend on rare events, whereas 
for large values of V, we expect that long simulations of large systems would be needed to 
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FIG. 7: (Color online) Critical rate w c versus parameter [i for T> = 0.02 and other parameters as in Fig. [6] 
Points: simulation; solid curves: diffusive MFT (upper) and simple MFT (lower). Error bars are smaller 
than the symbols. 



observe the crossover from mean-field-like behavior to the asymptotic scaling regime. 

The spreading exponents 9, § and z sp are estimated via analysis of the local slopes, 9(t), 
etc. For example, 9{t) is defined as the inclination of a least-square linear fit to the data for 
n{t) (on logarithmic scales), on the interval [t/a, at]. (The choice a represents a compromise 
between high resolution, for smaller a, and insensitivity to fluctuations, for larger values; we 
use a in the range 2-3.) We estimate the exponents by plotting the local slopes versus 1/t 
and extrapolating to 1/t = 0. Since a supercritical process leads to local slopes that curve 
upward, and vice-versa, we seek the value of w associated with the least curvature. The 
local slopes are plotted in Fig. [HJ On this basis of these results we estimate the critical point 
as w c = 0.2972(1), and find 9 = 0.573(5), S = 0.096(1) and z sp = 1.768(2). (The estimate 
for w c is based on the data for 9{t).) The results for the exponents compare reasonably well 



with the literature values, 9 
two dimensions 
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0.586, 5 = 0.092, and z sp = 1.771 for dynamic percolation in 



271 ]. The main source of uncertainty in the exponent estimates is the 



uncertainty in w c itself. The exponents obey the scaling relation of dynamic percolation, 



9 



•up 



25 — 1, to within uncertainty. 



To determine the fractal dimension of the critical cluster, we studied the radius of gyration 



R g of the final cluster as a function of its size, n, in a set of 500 realizations using tM = 5000. 
One expects that at the critical point, R g oc n l / Df , for n ^> 1 [25]. A least-squares linear 
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FIG. 8: Simulation: local slopes 6(t), 5(t), and z sp (t) for (lower to upper, on left) w = 0.297, 0.2972, and 
0.2975. Parameters: V = 0.02, A = 0.05, // = 0.5, v = 0.1, and wtr = wdr = wts = wds = 0.2. 



fit to the data (see Fig. [9]) yields R g ~ n°' 525 ^\ corresponding to a fractal dimension of 
Dp = 1.91(2). The value for two-dimensional percolation is 91/48 ~ 1.896 25]. 



B. Subcritical regime 



In the subcritical (nonspreading) regime, three quantities of interest are the mean lifetime 



t p , the mean (final) cluster size n and the mean radius of gyration R g of the fina 
One expects that, in the neighborhood the critical point, these quantities scale as 



cluster. 



t p oc A" 



n oc 



A" c , R s oc A" 



(21) 



where A = w c — w, and ( = 7^||/Vj_, with 7 the percolation critical exponent governing the 
divergence of the mean cluster size. 
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FIG. 9: Simulation: radius of gyration R g versus cluster size n at the critical point for parameters as in 
Fig. E with w = 0.2972. 



We estimate the exponents mi, ( and v± using simulations with system size L = 1000, and 
1000 realizations for each value of w studied, for the parameter set defined above (see Fig. 
HOl) . The simulation results yield the estimates u\\ = 1.52(2), u± = 1.29(3), and £ = 2.69(3). 
(We note, however, that the result for u± should not be taken very seriously, given the 
small values of R g .) The reference values for dynamic percolation in two dimensions are 



z/|l = 1.506, u ± = 4/3, and C = 2.698 |26|. 



C. Clusters and spreading 

Figure [11] shows examples of growing clusters at the critical point for two rather different 
values of the diffusion rate. The one corresponding to a smaller T> value is more densely 
connected, while the other shows evidence of "colonies" growing at some distance from the 
main concentration, as well as a more diffuse boundary. The distribution of transformed 
and depleted cells around the perimeter is highly nonuniform in both cases. Further growth 
appears to be likely only in limited 3X6 clS, £lS reflected in the nonuniform signal concentration. 
In the supercritical regime, cluster growth is more symmetric, but still somewhat irregular, 
as shown in Fig. [121 f° r parameters such that w ~ 1.08w c and 1.33w c . 

We close this section with results on the propagation velocity v s in the spreading phase. 
Near the critical point (or the critical line in the w-T> plane) the velocity is expected to scale 
as v s ~ e v w~ Vx , where e is the distance from criticality (sjj. This gives v s ~ e - 173 for dynamic 
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FIG. 10: Simulation: survival time t p , mean cluster size n, and mean radius of gyration R g versus A = 
w c — w, in the subcritical regime, for parameters as in Fig. [S] The slopes of the regression lines are 1.52 
(t p ), 2.68 (n), and 1.29 (R g ). 

percolation in two dimensions. In simulations, we determine the spreading velocity via the 
relation (Nn(t)) ~ ir(v s t) 2 , i.e., the region of removed cells tends, on average, to a circle of 
radius v s t at long times. In these studies we perform ~ 100 realizations for each w value, on 
lattices with L = 850 - 1000, extending to a maximum time of = 8000 units. The results 
(Fig. [TBI are consistent with a power law near the critical point. A fit to the data, using 
w c = 0.1206, yields v\\ — v±= 0.178(15); including the uncertainty in w c itself (±0.0001), we 
obtain v» — v^= 0.18(4), which, while not very precise, is consistent with the value expected 
for dynamic percolation. The inset of Fig. [13] confirms the scaling v s ~ VT>, as expected on 
the basis of dimensional analysis, away from the immediate vicinity of the critical point. 
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FIG. 11: Growing critical clusters at time 2000, for parameters as in Fig. [6l with T> = 0.02, w = w c = 0.2972 
(left), and T> = 0.3, w = w c = 0.0812 (right). Light color: removed cells; black dots: transformed cells; 
black diamond: position of original seed; x: positions of relatively high signal concentration, C > 0.1. 




FIG. 12: Growing clusters in the spreading regime. Contrasting colors show the set of removed cells at times 
100, 200, 500, and 1000. Parameters as in Fig. El with V = 0.1. (For these parameters w c = 0.1205(5).) 
Left panel: w = 0.13 ~ 1.08w c ; right: w = 0.16 ~ \.3Zw c . 



V. DISCUSSION 



We study an epidemic model consisting of elements (organisms in a community or cells 
in tissue) with fixed positions, in which disease or damage is transmitted by diffusing signals 
emitted by infected individuals. The model is formulated on a square lattice in which each 
site bears a cell, which can be in one of the states: susceptible, transformed, depleted, 
or removed. Signal diffusion and decay is treated deterministically, given the (random) 
source distribution in space and time. We study the model using mean-field theory (in both 
simple and diffusive versions) and Monte Carlo simulation. Simple MFT predicts the order 
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FIG. 13: Simulation: spreading velocity v s = (Nr) 1 / 2 /(y/nt) versus w for V = 0.1. The slope of the 
regression line is 0.178. Inset: spreading velocity versus V 1 / 2 for w = 0.15. Other parameters as in Fig. [SJ 
Lines are a guide for the eye. 



of magnitude of the critical value w c , if the diffusion rate is not extremely small, but is 
insensitive to the diffusion rate. Diffusive MFT yields a slight improvement over the simpler 
analysis; it captures, qualitatively, the fact that w c decreases with V, and that w c diverges 
as V ->■ 0. 

The process is found to exhibit a continuous phase transition between spreading and 
nonspreading phases. Simulations of the spread of activity yield estimates for the critical 
exponents 9, 5, and z sp consistent with those of two-dimensional dynamic percolation. The 
fractal dimension of the cluster of affected individuals at the critical point is also consistent 
with that of dynamic percolation, as are the critical exponents associated with the subcritical 
regime, and the scaling of the spreading velocity in the supercritical regime. Although these 
results are obtained for a specific set of parameters, there is little reason to expect a change 
in scaling behavior for other values, as long as the diffusion rate is finite. Indeed, dynamic 
percolation universality for epidemic-like processes with a finite range of spreading was 
already asserted some time ago sj]. We are unaware, however, of a previous verification of 
such behavior in the case of propagation via a diffusing, decaying signal. 

The present study suggests several lines for future work. One is a more detailed study of 
the scaling of the spreading velocity. The ability of mean-field theory or reaction-diffusion 
equations to describe this aspect of the process is of interest, as such approaches are fre- 
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quently used in applications. Another subject for future study concerns the nature of the 
spreading transition in disordered and fractal media. The possibility of a discontinuous 
transition for a nonlinear dependence of the transformation rate on signal concentration 
is also worth investigating. While mean-field theory does predict such a transition when 
the concentration-dependent rates are oc C 2 , experience with contact-process-like mode 



suggests that the nature of the transition depends on the details of the dynamics [28|, |29 |. 
Finally, applications to specific processes, such as the radiation- induced bystander effect, 
are of great interest, if plausible estimates of the governing rates can be obtained. 
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